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ABSTRACT 



The estimation of time delay and Doppler difference of a signal arriving at two 
physically separated sensors is investigated in this thesis. Usually, modified cross power 
spectrum coupled with Doppler compensation is used to detect a conunon, passive signal 
received at two separated sensors. Another successful approach uses the cross coherence 
to achieve this goal. This thesis modifies these two techniques to model the Doppler 
difference via an autoregressive (AR) technique. Analytical results are derived and ex- 
perimentally verified via a computer simulation. Performance at high and low signal to 
noise ratios ( SXRs ) is examined. 
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1. INTRODUCTION 



This thesis investigates the use of autoregressive (AR) models for estimating the 
Doppler difference and differential time delay by processing of a narrow band signal 
emitted from a moving source and received at two physically separated sensors. If the 
signal is received at two different geographical positions even in the presence of uncor- 
related noise, then, depending on the signal strength and duration, it is possible to esti- 
mate the differential time delay. 

Because the source is moving, the signals that are received at the sensors may have 
different frequencies due to the Doppler effect. To obtain accurate differential time de- 
lay estimation. Doppler difference compensation is usually required. 

This compensation can be implemented by using frequency shifting of the narrow 
band components of the received signal. This frequency shift can be obtained using a 
Fourier transform. In this thesis we use an AR model to detect the frequency shift. 
Using this Doppler compensation, an estimate of differential time delay can be obtained. 
Estimating the delay and Doppler using an AR model can be interpreted as a form of a 
narrow band coherence procedure, provided the estimations are properly normalized. 

This thesis is arranged in si.x chapters and six appendices. Because the estimation 
of the time delay and Doppler is intimately related to the coherence between two trans- 
formed complex signals, an extensive investigation of coherence is given in Chapter II. 
In Chapter 111. AR models, advantages of AR modeling, and AR model order selection 
are presented. Chapter IV is devoted to the analysis and the processing of noisy signals 
to estimate the differential time delay and the Doppler difference. To estimate the dif- 
ferential time delay, two approaches are pursued. AR model performance, Doppler es- 
timation and two types of time delay estimation are examined in Chapter V. In the last 
chapter some general conclusions of the work carried out in this thesis are presented, and 
some suggestions for future investigations are given. Computer simulation programs are 
included in Appendices E and F. 
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II. COHERENCE 



A. DEFINITION 

The coefTicient of coherence between two wide sense stationar}' random processes 
is the normalized cross power spectral density function defined by Wiener [Ref 1: p. 12] 
as 



Vxyif) = 



Gxyif) 

\ Gxxif)Gyy{J) 



( 2 . 1 ) 



where / denotes the frequency (Hz) , 

is the cross power spectrum between .v(/) and y(/) . 

G^XO denotes the auto power spectrum of x(r) . and 
G^X.f) denotes the auto power spectrum of r(t) . 

Despite some confusion in the literature. Wiener intended for the coefficient of co- 
herence to be complex. This is apparent since he discusses both the modulus and the 
argument of the coefficient of coherence. 

B. PROPERTY OF THE COHERENCE FUNCTION 

The power spectral density matrix 0(f) is positive semi definite. Therefore, for two 
random processes .v and y. we see that 



[0(/)] = 



GxxiO Gxyif 

GyxU) Gyyif) 



> 0 



( 2 . 2 ) 



For real processes we have G^,(f = GXU) nnd thus 

GxffGyyif)- \Gxy(f \^>0 



where * denotes the conjugate of a complex number. .And 

Gxxif)Gyyi/)> 



(2.3) 



(2.4) 



Furthermore, G,Jf) and Gyy(f are nonnegative, real functions of frequency. When 
G,,(/) and GJf) are strictly positive definite ( i.e., G,XfG^XJ) > ® ). Eq. (2.4) can be di- 
vided by G^XfGyXf) \vithout changing the sense of the inequality. 



This provides as an upper bound 



!>V(/)I<1 V/ (2.5) 

and since the magnitude of any complex number is greater than or equal to zero, we 
have the lower bound 



0< |y^,.(/)| <1 (2.6) 

The magnitude of the coherence function is always between zero and one. It is zero 
if the processes x{i) and y(r) are uncorrelated and it is equal to unity if there exists a 
linear relationship between x{r) and j(/) . In order to define the coherence it is necessaiy 
that the numerator and denominator of Eq. (2.1) are not simultaneously equal to zero. 
Coherence is not defined if either auto spectra is zero. 

C. COHERENCE ESTIMATION 

If .\’(7)) denotes the Fast Fourier Transform (FFT) of the / th segment of w(«) at 
freqtienc\\4 > then the spectral density estimates are given by [Ref 2: p. 22] 

V 

/=! 

Lm = ^ym)yi(fk) (2-s) 

Z_j 

/=1 




where a = ■ 



1 



XTf, ’ 

A = number of segments, 
T = segment length, and 
f, = sampling frequency. 



Finally the coherence estimation is 



(2.9) 
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m = 



\ 



^xyifk) 



xm'lifu) 



( 2 , 10 ) 



/=1 



y,i-vtt)i^yi w«)i' 



/=i 



i=\ 



D. COHERENCE OF NARROW BAND SIGNALS WITH DIFFERENTIAL TIME 
DELAY AND DIFFERENTIAL DOPPLER 

The output of tlte band pass filters (BPFl and BPF2) in Figure 1 are denoted by 
A'(/I) and }';(/^) respective!}'. Each term represents the Fourier transform of the corre- 
sponding time series evaluated at frequency and time /. For narrow band signals a 
Doppler shift corresponds to a frequency shift. If a signal arrives at the two sensors 
having a difierential Doppler shift as well as a dilTerential time delay, then we see that 
a frequency compensation and time delay compensation by the appropriate values tend 
to line up the signals in frequency and time. This is accomplished by using an additional 
Fourier transform in channel- 1 of the processor and a time delay in channel-2 of the 
processor. Mathematically, this can be expressed as 

V 

!'W) - (-Ml) 

V /=] /=i 

Comparing this with Eq. (2.10), we see that we generalized the coherence concept. We 
also note that the implementation resembles a correlator, where one of the signals is 
frequency compensated and the time delay corresponds to the delay operation of the 
correlator. 
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Figure 1. Coherence estimation block diagram. 



If the Figure 1 is redrawn as in Figure 2, then it can be interpreted as an FFT 
implementation as shown in Figure 3. 




Figure 2. Coherence estimation block diagram (reinterpretation). 
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Figure 3. Coherence estimation block diagram using the FFT. 



But this FFT lias a poor resolution. To obtain a better resolution, an AR model 
is desirable. This approach is shown in Figure 4. 




Figure 4 . Coherence estimation block diagram using an AR model. 



Figure 1 and Figure 4 represent two difFerent schemes to compute the coherence 
function. Dilferential time delay and Doppler dilFerence can be estimated by using AR 
models in the modeling of the coherence function. 
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III. AUTOREGRESSIVE (AR) MODELS 
A. AR MODELING 

If the following dilference equation is satisfied, the resulting structure is called an 
AR model of order p. [Ref 3: p. 177] 



p 

•^('0 = - -/<) + «(») (3.1) 

k=\ 

where x{n) = the signal at instant n, 
u{it) = the white noise driver, and 
o* = the /< th coefficient of the AR model. 

A realization of the AR model is illustrated in Figure 5. 




Figure 5. AR filter of order p. 



7 



B. ADVANTAGES OF THE AR MODEL 



The motivation for parametric modeling of random processes is the ability to obtain 
better spectral estimates based upon the model than estimates produced by classical 
spectral estimation. Both the periodogram and correlation methods can be used to yield 
Power Spectral Density (PSD) estimates. The unavailable data or autocorrelation se- 
quence (ACS) values outside a given window are assumed to be zero. This kind of an 
unrealistic assumption leads to distortions in the spectral estimate. 

The advantages of the AR approach are 

1. AR spectra tend to have sharp peaks, a desirable feature of high resolution spectral 
estimators. 

2. AR parameter can be obtained as solutions to linear equations. 

C. POWER SPECTRAL DENSITY 
Eq. (3.1) can be rewritten as follows 



p 




oo 




(3.2) 



k=Q 



where /i* = the causal filter impulse response. 
Let us take the Z -transform of Eq(3.2). 



X(z)^H{z)U{z) 



(3.3) 



Rewriting Eq. (3.1) as 



p 




(3.4) 



where Cq “ 1- taking the Z -transform gives 



Aiz)X{z) = L-( 2 ) 



(3.5) 



Eq. (3.3) and Eq. (3.5) can be solved to obtain H(z) 



S 



(3.6) 



H{z) 



I 

A{z) 



and hence 



A'(r) = -^ (3.7) 

A(.) 

The Z - transform of the output sequence {x{n)] is related to the Z - transformation 
of the input random process u{n) by [Ref 4: p. 56] 

6-(2)f*(4) 

p,^z) = xiz)X{ -V ) = = Puui=) rT“ 

A(z)A {-^) A(z)A (^) 

The AR power spectral density is obtained by substituting z = e’-'^ into Eq, (3.S) 
and scaling it by the interval T. 



^ARif) — T'Po 



\Aif)V 



= Tp. 



(J)aa e.(J) 



(3.9) 



where c. = f \ .e-'^zfT 

/■ 

<7 = C 1. a.J'’ 

p„ = variance of driving sequence. 

D. BURG'S ALGORITHM 

In practice, the autocorrelation is usually not available, so one must make an AR 
spectral estimation based on the available data samples. The simplest procedure to ob- 
tain an AR spectral estimate from data samples would be to produce estimates of the 
autocorrelation sequence from the data. These autocorrelation estimates would be used 
in lieu of the true autocorrelation sequence in the YULE - WALKER equations to yield 
the AR coefficients. However better results are obtained, particulary for short data 
segments, by algorithms that obtain the AR model parameters directly from the data, 
without explicitly estimating the autocorrelation function. 

The Levinson recursive solution to the YULE - WALKER equations relates the 
pth order AR parameters to the p — 1 th order parameters as given by [Ref 3: p. 21 1 ] 
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For /I = 1 to /2 = /> — 1. the reflection coefficients can be found by using the 
known autocorrelation function for lag 0 to — 1. So we have 

p-\ 

< 3 - 1 1 ) 

The recursion for the driving white noise variance is given by 

Pp = Pp.,{\-\kp\^) ( 3 . 12 ) 



where Po = r„(0) . 

But the .‘\CS is not available, hence we can not calculate the reflection coeOlcients. 
The Burg algorithm provides an estimate of the reflection coefficients which in turm are 
obtained through a least squares criterion. 

The fonvard linear prediction error is given by 



= -v(«) + 



^ c/p{m).\{n — m) 



m-\ 



(3.13) 



while the backward linear prediction error is given by 



4 (») 



r 

= x{n - c/p [m).x{n + m - p) 



(3.14) 



m=l 



Substitution of Eq. (3.10) into Eq. (3.13) and Eq. (3.14) yields the recursive relationships 



ep{n) = c^_,(/i) + kpcl_^{n - 1) 


(3.15) 


= Vi(« - 1) + kp(/p-M 


(3.16) 



At each order p, the arithmetic mean of the forward and backward linear prediction 
error power ( sample prediction error variance) is given by 
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(3.17) 




^ S Yj 



J(«) I ' 



This expression is minimized, subject to the recursion given by Eq. (3.15) and Eq. 
(3.16). Thus, pf is a function of single parameter, namely the complex valued reflection 
coefficient k^. Setting the complex derivative of Eq. (3.17) to zero 



Spb , . gpf 
cRe{k^) • J dlm{kp) 



(3.1S) 



and solving for k. yields 



A 




-2 V - 1) 

n^p-k~\ 



(3.19) 



The estimation of the reflection coelTicient represents the HARMONIC mean be- 
tween the forward and backward partial correlation coefficient, where 
e^{n) = 4(i;) = .v(/?) and .V == number of data points. 

E. FINAL PREDICTION ERROR (FPE) CRITERION 

Because the best choice of filter order is not generally known a priori, it is usually 
necessary in practice to postulate several model orders. FPE is a kind of criterion which 
was provided by .Akaike [Ref 3: p. 230]. This criterion selects the order of the AR 
process so that the average error variance equals the sum of the power in the unpre- 
dictable part of the process and of a quantity representing the inaccuracies in estimating 
the AR parameters. The FPE for an AR process is defined as 



FPEip) = p. 



N+(P+1) 

A-(p+l) 



(3.20) 
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where A’ is the number of data samples, 
p is order of the filter, and 

Pp is the estimated white noise variance when using a p th order filter. 
The order p selected is the one for which the FPE is minimum. 
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IV. DOPPLER AND DIFFERENTIAL TIME DELAY ESTIMATION 



A. DOPPLER ESTIMATION 



Let x(/) and ;•(/) be the signals received by two sensors 



■v(t) = cos{2;r(/+ a,)/} 



(4.1) 



y(r) = cos{2n{f + aj)/} 



(4.2) 



where /is the carrier frequency, 

a, . a, are Doppler shifts, and 
i is the time variable. 

Let / be the sampling frequency which satisfies the Nyquist theorem. For con- 
venience. in all derivations and simulations a sampling frequency of 64 Hz is used, to- 
gether with band pass filter width of 1 Hz. We assume that |a, | <0.5 ( /= 1,2 ) . 
which implies that the signals stay in the band pass filter regions of their respecti^'e band 
pass filters regardless of any Doppler shift. Note, these values can be modified to arbi- 
trary sampling rates and pass band regions. The sampling rate is given by 



/>2/+ l>2(/-f ^,) ( /=1.2 ) 



(4.3) 



The phase at instant k and sampling time interval T are given by 




(4.4) 



^ 2n{f+ &. 2 )k 

y4>k - 7 



(4.5) 




(4.6) 



Let x{n) denote the sampled analog signal .v(t) | then 



.y (/2 -f m) = cos{2n(f + a,) ^ -f 

Js 



(4.7) 
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(4.S) 



y{n + m) = cos{2n{f + c^..) ^ 

Js 

The derivations ofEq. (4.7) and Eq. (4.8) are given in Appendix A. 

These signals. (.v(/7)} and {.v(«)}. are processed to detect the Doppler difTerence . 
Let BPFl be the band pass filter centered at J] and BPF2 be the band pass filter centered 
at/, (where / might equal/ ). A’ data points are processed in the band pass filters, to 
generate one output. The inputs to BPFl and BPF2 at time m are vectors A'„(/?) and 
}'„(//) . respectively. The size of the input vector to the filter is the number of data points 
taken during 1 second (i.e.. the number of points processed in the filter = .V=/ ). The 



input vectors are denoted by 

A/(«) = [-v(m). -v(//; + 1) x{m + A' - 1)] (4.9) 

•;■(/« + 1) y(//i + .V — 1)] (4.10) 



The filtering is performed using FFTs . where successive EFT outputs are generated at 
the input data rate. The BPF'2 output is conjugated. 

To avoid the complexity the following four complex constant variables are defined. 

v-i 

.4^= ycos{2;rf/+a,)y (4.11) 



.V-1 




.V-1 

A, = V cos{2;t(A+ «2) f (4.13) 

r.=0 

v -1 

By = sin(2;:(/' + « 2 ) v (4. 14) 

At instant ni the complex output signal of BPFl can be calculated as follows 
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V-1 



= ^cos(2;r(/' + a,) ^ 

;;=0 



V -1 



=I 

/ 7=0 



[ cos{27r(/ + a,) } cos - sin{2;r(/' + a,) ^ } sin Jc a' 

Js Js 



A -1 A -1 

= coSj,(/)„,^cos{2r(/ + sinj,(jf)^^sin{2:r(/' + a,) 

/,=o " «=o ^ 

= cos sin 

_j_ ^ J m ^ m m 



= ''^.X -); 



+/-^x J,<?„, , ^-x 2 ^x _-, 

^ t “T ^ 






4-//J .^ . {/+ =<i) /-/ —iB ■ C^+*i) 

^.v jlo — T m J^\x -j2(p m 

^ A ' ^ Js 



The complex signal XJf) contains two frequencies with different amplitudes. 

To understand the character of A'„(/) , it is important to evaluate which is the 
dominant term in the above expression. 



v -1 



+J^X = 7 c cos{2;r(/' + a,) — } +j sin{2;r(/’ + a,) — -V 

^ Js Js 

n=0 
V -1 






/i=0 

v -1 






«=o 



1 - 

1 



( 4 . 16 ) 
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.V-] 



A -J^x = cos{27t(/ +ax)jr]-J sln{2z(f + a,) jr }]^ ^^ 

n=0 ^ 

A’-l 






/!=0 

.-V-1 



- Jl -{ 2f + x0jr 



(4.17) 



n=0 



1 - 

1 _ 



From Appendix C 

\A^+jB^\>\A^-jB_^\ (4. IS) 

Therefore is the dominant term of .V^(/). 

In simiiiar way, the complex output signal of BPP^2 can be calculated as follows 



v-i 

}•;(/) = ^cos{2;r(r+ a.) f + 

K^-0 

.V-1 

= cos{2nlf + a,) ^ } cos y<t>„, - sm{2n(f + 0 . 2 )^} 

/:=U 

.V-1 \'-l 

= cos v,0„,T cos{2;:(/'+ a.j) -j- }e^'~ ~ - sin v</)„,^sin{2r:(F + ^ 2 ) A 

(4.19) 

= Ay COSy < f )„, - By SIH 



= A .. 



^ jy <? rr , £ J 






Ay +/^v • • ~~jBv • • 

= -J— ..-j, >. 

^ c r 2 ^ 

A >+ J^y A ~ J^y - n ^ 






e - n < P—^rn 



The complex signal Y'JJ) contains two frequencies with different amplitudes. 
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Even though Y’„{f) is siiniliar in form to XJJ), it is not obvious which term is dom- 
inant. 



A-l 



Ay +jBy = cos{2;:(/'-i- ^ } +j sin{2;r(/' -t- aj) ^ 'V 

rt =0 ^ ^ 

A'-l 



/?=0 

A-l 

j _ ^/2:r(2/+ a,) 

1 _ ^/2:r(2/+:.,)-^ 

, J2t5, 



(4.20) 



1 _ ^- 2 ^( 2 /+ 



A 



A-l 



JL mJM-t 



y J^y- . [cos(2-(/-+a,)-^}-ysin{2-(/'+a,)^}]/ 

Js Js 



n={) 

A-l 






n=0 



A-l 




«=0 



- 22-32 



1 - e A- 



(4.21) 



From Appendix C 



I Ay —jBy \ > \ Ay +J B y \ 



Therefore 




is the dominant term of Y'JJ). 



(4.22) 
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Since the two output sequences {A'„(/)} and { ¥'„{/)] have two sinusoidal components 
each, their product {XJf)Y'^.{f}} has four sinusoidal components. The four frequencies 
are 



cui = 2n 



(Ki - O',) 



(u, = 2;r 



(y-2 - °^i) 

fs 



, (03 = 



(2/ + g| +0:2) 
fs 



, W4 






ITT- 



(2/+«i +0:2) 



f 



with cui the sinusoid with the largest amplitude and a >2 the sinusoid with the smallest 
amplitude. 

The product of the output of the filters is given by 



a;, (/)>;*,(/) 



'^x+J^x ''f J^y e 



{Ax -rjB^){A^. —JB^) _ 

' pl^x-m m> 



+T(ai. y-2. x4>m-y<l^m) 
+T(a,. y- 2 - x*^m- y4>m) 



(4.23) 



where F[y.^, y.^. .4>„. ,</>„) represents the three low amplitude frequency terms. 

\Mien using the AR model as described in chapter III, at a high SXR ( i.e., 
S.\R-*oo ) . a 4th order AR model detects all of the four frequencies given above. 
When the S.V/? is low (i.e.. 20 dB ). only one dominant sinusoid is detected with fre- 
quency 



fa), = • 



IrAf^y.^) IrAf+y.,) 



f f 

(ai - a,) 

~~7s 



(4.24) 



Since we can detect co, and know the value ofy^ , a, — aj can be estimated by using an 
,AR model. 

B. DIFFERENTIAL TIME DELAY ESTIMATION 

To detect and localize a target, it is important to find the diflerential time delay of 
signals arriving at two sensors. Let us assume that signals at the two sensors are as 
given in Figure 6 ( i.e.. zero dilTerential time delay ) For the remainder of the figures, the 
time axis is scaled to be 64 points per second. 
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Figure 6. Receiving signals at sensor 1.2 (SNR= 100, no differential delay). 



We can estimate the difl'erential time delay and differential Doppler using two dif- 
ferent approaches. The first approach uses the cross power spectrum while the second 
approach uses the coherence. 

1. Differential time delay and differential Doppler estimation using the cross power 
spectrum. 

Figure 4 shows the block diagram of the coherence estimator using AR mod- 
eling. In this figure, the output of the FFT can be interpreted as the cross power spec- 
trum. Using this cross power spectrum, we can estimate the differential time delay and 
the differential Doppler. But when we use the highest peak of the cross power spectrum, 
the peak is somtimes not detected at the proper time delay nor at the proper Doppler 
shift. We will show a special case in which the peak of the power spectrum is not de- 
tected at the proper time delay nor at the proper Doppler values. 
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a. Special case. 

For our test case the data duration is six seconds, two seconds of the noisy 
signal and four seconds of noise only. Linear transformation of this data leads to one 
of three types of outputs. The first type represents full information, while the second 
type represents partial information. The third type represents a noise only condition. 
Generally when two signals are lined up in time, the .A.R model should give the highest 
output power. But this is wrong in some cases. In the full information case, the mag- 
nitude of transformed signal is high and constant. For some reduced information case, 
even though the magnitude of the transformed signal decreases, it still may have large 
magnitude of spectral components. This phenomenon can be explained as follows. 

Define two functions depending on k 



S-\-k 



j,Ay= cos(2t:/2 



2r/^ 



n=0 



Vs 



V = y sin(2-/2 y A- 



(4.25) 



(4.26) 



where =/+ oc;,. 

k denotes the number of lost data at BPF2 input vector. 



We assume that we know when the BPFI signal starts. If this information 
is not available, we need to exanrine the signal at the output of BPFI. to obtain a can- 
didate time frame. Consider the input of AR model as shown in Figure 7. 

The input data size to the AR model is the number of linearly transformed 
data points during a given period (i.e., input size is A= f , ). The BPFI output represents 
full information (i.e.. each output element is produced from the signal information 
points). Therefore, as shown in the previous section, the dominant output sequence is 

A-,y. kV‘^\ , kV'^-'-'-' (4.27) 



where 




1 - 



(4.28) 
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Figure 7. AR model and its driver source. 

f ] ) 0 second delay case. We already defined tlie input vector and 

)'„(/;) to BPFl and BPF2. This vector can be interpreted as a time snap shot. Fach 
output of the BPF requires /V input points. Ifwe require A' output points from the BPF, 
even using maximum overlap in the processing, we require at least 2A'— 1 points at the 
BPF input. Figure 8 shows the 2N input to BPFl and BPF2. Both IN input points 
contain tlie signal and some noise. As shown in Figure 8, in the zero delay case, the two 
input vectors X„{n) and F„,(/?) ( 0 < < N ) provide full information. So the lin- 

early transformed output XJJ\) and Y'Jf^) also represent full information. 

The dominant part of the BPF2 output is the sequence 

(4.29) 



where 






1 1 - 






(4.30) 
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Figure 8. T^^o signals. 0 second delay at the sensors. 



The dominant input sequence to the AR model is given by 

.A- = 0,1.2,...,.V- 1 } (4.31) 

where the magnitude is A, A, 

flj -/ second delay case. Figure 9 shows the 2.V input to BPFl and 
13PF2. The 2.V input data points of BPFl contain the signal and noise. During the first 
.V points the input to BPF2 contains the signal plus noise while during the second .V 
points only noise is present. The input vector A'„(«) of BPFl contains full information, 
but input vector of BPF2 does not contain full information. In the }„(a 7) case, 
every element of To(/0 full information. But when A is not equal to zero then >*(/;) 
consists of .V— A information elements and A noise elements. 




Figure 9. Two signals. -1 second delay at tlie sensors. 



The output from BPF2 represents partial information while the out- 
put from BPFl represents full information. The BPF2 output sequence is given by 



{ 



A-4 +JkBy ^ k^y JkBy 



k = 0A....X- 1 } 



(4.32) 



The derivation of Eq. (4.32) is given in Appendix B. 



Assumins that 



l,Ayt jBy 



^ is the dominant term, which is a valid assumption for our 

test case with j^ = 64//z, f=23Hz, a, = 0.4/Fz and aj = 0.001 //z, then the input se- 
quence to the AR model is dominated by 



{/c, k = 0,l.2, A'-l} 



(4.33) 



with ki is defined earlier ( Eq. (4.2S) ). 
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The magnitude of dominant input sequence is obtained by examining 



,V-I-Ar 



y [cos{2-{f+ a,)} -J sm{2n(f+ 



n=0 

\-\-k 






n=0 

X-\-k 



(4.34) 



-P-^2 — 



;?=0 

n 

1 - 



When we compare the magnitude of the sequence in Eq. (4.27) and the magnitude of the 
dominant term in Eq. (4.32). the magnitude of l/tj > I ^.4, — | /2 where 



A-. = 



I 1 -c 



—J2r:OL2 



- 1 _ ,v 



and 



\'-k 

Pp'^y ~Jk^y 1 1 — ^*2 



1 _ ^ y 



.40 



,-/2^S2 



1 — e 



—J 2 -X 2 



-L 2 






1 — e 



-72 - 3 <; 



(4.35) 



\\'c note that Eq. (4.35) has two frequencies ( i.e.. 0 and 2z ). 



The masnitudes of both terms are 



1 - 



Now. note that as long as |aj| < . then the magnitudes of both terms in Eq. (4.35) 

are larger than I Ay I 



Due to the symmetr}' a positive delay results in a similar reponse to a negative delay, and 
hence Eq. (4.35) holds also for delay of + 1 second. Equation (4.35) can be interpreted 
as having large responses at shifted Doppler values at a delay of 1 second and a delay 
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of — 1 second. Proper delay can not be established nor can the Doppler value be estab- 
lished. This cross power spectrum algorithm will not give good information about either 
Doppler diflerence or time delay. 

b. Modified cross power spectrum. 

To detect the signal, input signals to the AR model should be preprocessed 
as shown Figure 10. If the BPF2 output signal is magnitude normalized, then this nor- 
malized signal retains good phase information. When we normalize the BPF2 output, 
one of two conditions can occur. The first condition ( i.e., synchronized ) leads to the 
magnitude normalization of Eq. (4.29) and provides accurate frequency and delay esti- 
mation. The second condition occurrs when the signal in channel-2 is not lined up in 
time with the signal in channel- 1 ( i.e., during the delay search ). Under this condition 
smaller peaks in the time Doppler plane appear at incorrect values of time delay and 
Doppler. Above 70 dB S.VR, the correct peak is the dontinant one and allows proper 
estimation, information from contour plots can be used at values of SA'/^ between 70 
and 10 dB ( see Appendix D ). 




Figure 10. Modified cross power spectrum block diagram. 



2. Differential time delay and differential Doppler estimation using the coherence. 
Another way detecting the signal is normalization of the AR output with the 
squared root of the product of P„(/) ( auto spectrum of x ) and P^Jf) (auto spectrum of 
y ). This can be done using two additional AR models as shown in Figure 1 1. 
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Figure 11. AR model of coherence. 



Using the AR approach, p„ , , p^, , and AR coenicienls can be ob- 

tained. Therefore we get 



Px,(/) = 



Tpxx 

\Axx^\^ 



(4.36) 



PyyiJ) = 



TPyy 



Ay, if) I ■ 



TPxyxy 

‘ xyx,\J) 2 

i A^y_^yij) I 



(4.37) 



(4.38) 
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(4.37) 



Pyyif) = 



TPyy 

I Ayyif) I ' 



Pxyxyif) = 



^Pxyxy 

M.,xv(/)l' 



where P„(/) is the power spectrum of A'/ , 
is the power spectrum of , 

P.yxJJ) is the power spectrum of , 

is the AR power spectral density of A', , 

Ayyif) is the AR power spectral density of . 

AryxyW is the AR power spectral density of A',}';^^ cross term. 

p,, is the driving noise variance of A’/ . 

pj.y is the driving noise variance of Y /.„- . 

p,„,, is the drivins noise variance of A; 17 _^ . and 

J S 



(4.38) 



We obtain a coherence estimate from the three power spectra, by assuming that 



\P.y{tW-\PxyxyijM • 



l^.x.xl/)/7vv(/)l 

\Pxyxyi/)\ 1 Pxvxv I ■4.x.x(0-4vvC0 I ^ 

"" \Pxx(f)Pyyif)\ ~ T l•4xy.xv(/)l 



(4.39) 



P 

1? 




Pxyxy 

PxxPyy 



M..x(/)4.vWl' 

Uxyx>(/)l 



(4.40) 



This approach provides good information about the differential time delay and 
differential Doppler down to SXRs of 20 dB. 
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RESULTS 



This chapter presents graphical results obtained by applying the analytical results 
of the previous chapters to specific examples and carrying out the required computations 
on the computer. 

A. AR MODEL 

This section provides the AR model performance tests. As we discussed in chapter 
111, Burg's algorithm is applied to this AR model. The FPE criterion is use to find the 
AR model order. Figure 12 is the power spectrum of two test signals sin(co,7'«) and 
cos(fc),r«) . Figure 13 shows the power spectrum of two other test signals sin(tU 2 r«) and 
cos(o),r/7). Both sets of test signals in these figures have an S.\R of 20 dB. The fre- 
quency/ is 13.45 Hz and/ is 23.45 Hz and T is second. 

Theoretically, the power spectrum of these signals have two impulse functions at/ 
and — / ( / = 1 ,2 ). Figures 12 and 13 show experimental results agreeing with the 
theoretical results. Figures 12 and 13 indicate which frequency has high power but they 
provide no phase information. 

ijrhe AR model is sensitive to frequency but is not sensitive to the phase. Fo r anv 
selection of frequenc}’ f which satisfies the Nyquist theorem, the AR model provides 
good frequency estimation provided the SXR is suflicient large. 

t 



:s 



INPUT SIN(cjTN) F=13.45 T=1/64 




wnai03ds d3M0d 



o 

CM 



M 

X 

o 

o 2! 
LU 
Z) 

o 

LU 

q: 



a 

CM 

I 



CD 



LO 

ri 



II 

ii_ 



tn 

o 

o 



X) 

CL 



1_ 

0009 




- o 




0 

- CM 

1 



I I I I 

000> OOOZ 0 

Wfldi03dS d3A\0d 



Figure 12. AR model performance test 1 (SNR- 20 dB). 
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Figure 13. AR model performance test 2 (SNR = 20 dD). 
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B. DOPPLER ESTIMATION 

In section A of chapter IV, we discussed five items. 

1. BPFl and BPF2 outputs have two sinusoidal components each. 

2. For BPFl. the donrinant frequency is designated as co, and for BPF2, the dominant 
frequency is designated as to, where 







-27z(f + a,) 



3. The power spectrum of cross term of BPFl and BPF2 has four components which 
can be detected in very high SXR (i.e. S.VR = 100 dB). 

4. When the S.\'R is low (i.e. SXR = 20 dB). only one dominant frequency compo- 
nent is detected. 

5. The dominant frequency of the power spectrum P,y,y{f) corresponds to the differ- 
ential Doppler frequency. 

Figure 14 shows the power spectrum of the BPFl output when /is 13 Hz. a, is 
0.45 Hz, and / is 64 Hz. Figure 15 shows the power spectrum of the BPF2 output 
when /is 13 Hz. a, is —0.45 Hz. and / is 64 Hz. As expected both figures show two 
spectral lines when the S.\R is 100 dB. The dominant frequencies are located at / equal 
13.45 Hz for BPFl, and at/ equal —12.55 Hz for BPF2. Figure 16 shows the power 
spectrum of the cross term of BPFl and BPF2. This figure shows four spectral lines at 
an SXR of 100 dB. while only one dontinant frequency is detected at an SXR of 20 dB. 
Dominant frequencies are always located between —1 Hz and 1 Hz. Figure 17 shows 
a subplot of Figure 16 for frequency between —1 Hz and 1 Hz. This figure shows the 
dominant frequency /= .9 Hz and the smaller spectral component at /^ —.9 Hz for an 
SXR of 100 dB. Only one dominant frequency can be detected at /equal 0.9 Hz for an 
SXR of 20 dB. Doppler difference from the dominant spectral component corresponds 
to the true Doppler difference a.i — exj which is 0.9 Hz 

Figure IS shows a comparison for SXR of 20 dB and 10 dB. Figure 19 is subplot 
of Figure IS for —1 Hz <f < \ Hz. From this figure, when SXR = 10 dB, the dominant 
frequency is shifted a little bit from the true Doppler difference. We have shown the five 
issues as we have discussed in the previous chapter. In summarx’, for S.VRs greater than 
or equal to 20 dB. the power spectrum of the cross term provides good Doppler esti- 
mation. 
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Figure 14. Po^^er spectrum of the BPFI output. 
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Figure 15. Pouer spectruin of the BPF2 output. 
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Figure 16. Cross po^^er spectrum of BPFl and BFT2. 
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Figure 17. Subplot of the cross pouer spectrum. 



35 



FREQUENCYfHZl 



SNR=20 





gOlxQ gOlXfr gO I Xf gO I X2 gO I X t 0 

wnHi03ds y3/.\od 



Figure 18. . Comparison the Doppler estimation at tno different SNR's. 
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Figure 19. Subplot of the Doppler estimation. 
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C DIFFERENTIAL TIME DELAY ESTIMATION 

In Chapter IV section B.l.a we discussed the problems associated with using un- 
normalized quantities and suggested two possible approaches to solve the problem. In 
this section we will show that for the case of ia, | <-^ Hz , the power spectrum of a 
nondelayed signal is smaller than the power spectrum of a signal delayed by one second 
{ i.e.. see discussion following Eq. (4.35) ). 

To estimate the time delay, we suggest the two different approaches. The first ap- 
proach is the normalization of the BPF2 output. The second approach normalizes the 
cross term P,y,y{f) with the power spectrum of BPFl T„(/) and the power spectrum of 
BPF2 Pyy{f) through the equation defined by Eq. (4.40). 

The following figures are simulations for /= 23 Hz, «i = 0.23 Hz, = —0.02 Hz. 
T = -^ second, delay ofO and an SSR of 100 dB. Figure 20 shows the surface plot of 
P,„SJ) Miile Figure 21 shows the contour plot of P,y.,y(J). Figure 22 shows the cross 
section plot of Figure 20 with respect to the DEL.AY axis. This figure shows that the 
power spectrum has a spurious peak at delay of one second. .-\t the proper delay (i.e.. 
0 second delay) has a relatively small value. 

P.vmXO is affected b}' two terms, namely and -r ~, — ^-rrr • Figure 23 shows the 



I J./'lP 

maximum power spectrum term of the transfer function and figure 24 shows the driving 
noise variance. At a delay of ± 1 second and of 0 second the transfer function shows a 
peak while the variance of driving noise has small value. In this case. P.y^JJ'), the product 



of the form 



at 0 second delav has a smaller value than that at anv other 



l-W/lP 

delay location. Using the scheme shown in Figure 10. we normalize the BPF2 output 
and present the result in Figure 25. ^^'e can detect the proper time delay and the proper 
Doppler difference. Figures 26 and 27 show the corresponding contour plot and the 
power spectrum of the normalized cross term. For any value of time delay and Doppler 
diflerencc. this approach provides good information about differential Doppler and dif- 
ferential time delay provided the SHR is greater than or equal to 70 dB. 

Results of the second approach { Figure 1 1 ) are demonstrated in Figures 28 and 
29. This approach also gives good information about differential time delay and differ- 
ential Doppler. This approach provides good time delay and Doppler estimation for 
SSR greater than or equal to 20 dB. 

If we use the contour shape, we can also estimate the differential time delay and 
diflerential Doppler for SXR of less than 20 dB. Differential time delay and differential 
Doppler estimation can be extended down to S\'Rs of about 0 dB when using the con- 
tour of the coherence surface. Figure 30 is a contour plot of using this second approach 
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( i.e.. coherence approach ) for a, = 0.23 Hz. «2 = — 0.02 Hz, Delay = 0, and an S.\R of 
0 dB. The contour shape of the second approach has the shape of the letter X. When 
the SXR is greater or equal to 0 dB, the proper time delay and Doppler difference seems 
to be located at the cross over point of the X. When the SXR is less than 0 dB, the 
contour shape does not shows an X clearly. In the coherence approach. Figure 31 
provides the contour plot, with channel- 1 containing only noise and channel-2 contain- 
ing signal plus noise ( SXR = 0 dB ). Figure 32 shows the contour plot, when channel-2 
contains only noise and channel- 1 contains signal plus noise ( SXR = 0 dB ). 
Figure 33 shows the contour plot, when both channels contain noise only. The above 
three figures do not show an X clearly, hence do not allow signal detection for esti- 
mation of any parameter. Using this information, we can distinguish the signal combi- 
nations from the noise combinations. 

\\'hen we estimate the time delay and Doppler difference using the contour plot, the 
modified cross power spectral approach has some problems as discussed in 
.Appendix D. 
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Figuie 20. Surface plot of the power spectrum. 
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CONTOUR PLOT OF POWER SPECTRUM 




Figure 21. Contour plot of the power spectrum. 
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Figure 22. Maxiinuni power spectrum. 
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Figure 23. Maximum power spectrum of the transfer function. 
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Figure 24. Variance of the driving noise. 
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Figure 25. Maximum modified cross power spectrum of the transfer function. 
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Figure 26. Power spectrum of the transfer function (surface). 
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Figure 27. Power spectrum of the transfer function (contour). 
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Figure 28. Surface plot of the coherence. 
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Figure 29. Contour plot of the coherence. 
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Figuie 30. 



Estimation using the contour shape (lou SNR). 
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Figure 31. Contour plot (noise in cliaiinel-1 only). 
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Figure 32. Contour plot (noise in channel-2 only). 
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Figure 33. Contour plot (noise only). 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



In this thesis, the din'erential time delay and dilTerential Doppler are estimated using 
two different approaches. The first approach uses a modified cross power spectrum al- 
gorithm while the second approach uses a coherence algorithm. In both cases 
autoregressive modeling of the required cross spectral component is used. The difi'eren- 
tial time delay and Doppler difference can be estimated using a threshold at high SXRs 
and a contour plot at low SXRs. 

At high SXRs ( i.e., S.\R > 70 dB ). the first approach locates the dominant peak 
at the proper time delay and the proper Doppler. The second approach utilizes two 
additional .AR models and two additional FFTs to obtain the autopower spectrum of 
channel- 1 and channel-2. At moderate SXRs ( i.e., SXR > 20 dB ), the second approach 
has the highest coherence peak at the proper time delay and Doppler frequency. When 
the SXR > 0 dB. the contour plot of the coherence has the shape of the letter X. The 
dilTerential time delay and Doppler difference can be located at the cross over point of 
the X. In high S.\Rs we can use the highest peak to detect the time delay and Doppler 
frequency using either approach. When we use the contour plot, the information about 
time delay and Doppler depends on the location and the form of the cross over point 
of the X. 

We can estimate the differential time delay and differential Doppler using the AR 
modeling of coherence. But we can not get the numerically correct value of the coher- 
ence coefficient. The following three points are recommended for future study. 

1. To get a cross power spectrum, we assumed that I I ~ I R,,!/) I* . This was 
required since we can not get a cross power spectrum using an .-\R model. Im- 
provement of the cross power spectrum estimates using an AR model can lead to 
improvement of the algorithm. 

2. Our first approach is not normalized by the auto power spectra. If we find the 
corresponding auto power spectra, we can also obtain a coherence coefficient esti- 
mate. 

3. To get a better information at a low SXR, we need to study the characteristics of 
the contour plots. 



54 



APPENDIX A. PHASE DERIVATION 



Define the phase at instant k. 



K^k - ' 



fs 



y4>k - 



2r.(f+ y. 2 )k 
fs 



Let x{n) denote the sampled analog signal x(/) | 
then 



.y(«) = cos{2n{f + K]) -^ } 

Js 

= cos{2n{f + a,) ^ + x4>q} 



x{n+ 1) - cos(2-(/ + a,) } 

Js 

fj H“ ^-i) 

= cos{2;:(/-+«i) J + ' ) 

= coi{2n{f + aj) -^ + y/)]} 

Js 



x{n + 2) = cos(2-t(/ + a,) ^ " } 

Js 

„ 2-{f+a,)2 

= cos{2rT(/'+ a,) f + } 

Js Js 

= cos{2r.(f+y.^)-^ + ^4>2} 

Js 



x{n + m) = cos{2-(/ + k,) ^ } 

Js 

,7 2nif-^a^)m 

= cos{2n{f + a,) — + j } 

Js Js 

= cos{2;r(/'+ + 



Similiarilv 



(4.4) 

(4.5) 

iciA) 

ia.2) 

(a. 3) 

(4.7) 



55 



r(«) = cos{ 2 ;t(/' + y.2) -y } 
Js 



= cos{2,t(/+ a^) ^ + >.'*0} 
Js 



y{n + 1) = cos{2;:(/' + a. 2 ) ^ ^ } 

Js 

^-2) 

= cos{2n(f + a 2 ) -^ + 7 } 

Js Js 

= cos{2;r(/'+ a 2 )Y + y<i^ 2 ) 



y{n + 2) = cos{2rr(/' + a 2 ) t " } 

Js 

n 2-(/'+a02 ^ 
= cos{2;rC/’+ a 2 ) ^ ~ } 

= cos(2;r(/+ 

Js 



y{n + m) = cos{2-(/' + a2) ~ } 

. n , 2-(f+y.2)m 

= cos{ 2 -(/ + ao) -y + 7 

Js Js 



= cos(2-(/' + a.2) -y + y 4 >f„} 
Js 



{a.5) 



ia.6) 



[a.l) 



(4.S) 
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APPENDIX B. OUTPUT OF BPF2 



Define two functions depending on k 

A'— 1 —k 

j,Ay= cos(2t:/2 



■S-\-k 

k^y= Y sin(27r/2^)/ 
^=0 



n.ff 



(4.25) 



(4.26) 



where = /+ a,. 

A. FULL INFORMATION CASE ( N POINTS ) 
.V points are data points at BPF2 input. 



,v-i 



yl,i/) = y cos(2;t/2 f + 

, 7=0 



V-1 



= cos{2nf2 jr ) C0Sy4>„, - sin(27T/2 y ) sin^0^}e'"’‘^ A' 



/2-/^ 



?i=0 



fs 



= COS 



A-1 A-1 

^ - sin^,(/)„,y sin(2rt/2 

n=0 



/7=0 



q / 1 ^. cos y(j) ^ O^y m 









~ O^y 






+7*0 






2y 
4 



(^•1) 
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B. PARTIAL INFORMATION CASE ( N-1 POINTS ) 
.V — 1 points arc data points at BPF2 input. 



.V-l-l 

}•;(/)= ^ cos(2;:/, + 



n=0 

.V-l-l 






{ cosilTzfj -Y ) C 0 Sy 4 )^ - sin(2rr/2 ) sm \ 

Js 



n=0 



fs 



.V-l-i 



.V-l-l 



= cos y4>„, cos(2r/, .V -siny4>^ ^ sin(27:/2y)c .v 



n^Q 



n=0 



fs 



(b.2) 



= iT. cos^,(p„, -,i?, sin 



m 



-vT 2 2J 

1 ^ V +J\ K i ^ 1 V ~J\ A' _ / 



, yy<?r- ^ 






C. PARTIAL INFORMATION CASE ( N-K POINTS ) 
.V— k points are data points at BPF2 input. 



X-\-k 



y„{J)= X cos(2;r/2^ + ,c/,„y-^'.v 



n=0 

\-l-A 



= { cos(2ry; -Y ) cos - sin(2;r/2 ) sin .v 



/ 7=0 



y: 



f 



\-\-k 



X-\-k 



= cos^,(j^„, Y cos( 2-/2 V — sin ,.(/)„ sin(27r/2 )c v 



;:=0 






«=0 



f 



(b. , 



= kf cos,(/),„-*B, sin (/) 



yr' m 



— j^A 















k^'^v . k‘f~Jk^v _/M 

— ^ ^ e ™ 



58 






APPENDIX C. AMPLITUDE COMPARISON 



Let C, and Cj be the complex amplitude. 






1 - 



Q = ' 



1 — e 



—jlnoi 



l - ,-^2.(2/+.)^ 



(C.l) 



(C.2) 



where | a 1 < 0.5, and N>2f + \ > 2{f+ a). The numerators of | C, | and | I are same, 
because the magnitude of the complex conjugate is the same. If we show that 
1 1 — I < 1 1 — , then we can say 1 C, | > 1 Cj . As shown in Figure 34 

1 1 - e-2-il = 1^1/1 and 1 1 - = \BI\ . 




The smallest angle LlOB is 
For positive f 



a) 

iV 



or 



2n(N-2f- a) 
iV 



2;r(l — a) 
N 



L10B = 



2rt(2/'+ a) 




2710. 


N 


> 


N 



LIOA 



(c.3) 
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Since 1 1 - a I > U | . 



LIOB = 


2n{ 1 — a) 


> 






A' 




X 



= L10A 



From Eq. (c.3) and Eq. (c.4) \AI\ <\BI\ , therefore 

I C, I > I C, I 



(c.4) 



(c.5) 



60 



APPENDIX D. CONTOUR PLOTS OF MODIFIED CROSS POWER 

SPECTRUM 



For SXRs greater than or equal to 10 dB the contour plot of the modified cross 
power spectrum algorithm has the shape of a boomerang. The shape of the boomerang 
depends on oc,. Figure 35 and Figure 36 show the contour plot of the modified cross 
power spectrum. The two figures show the boomerang shape and their different di- 
rections. For positive Kj ( case I ) and negative a, ( case 2 ) the opening of the 
boomerang is toward the left and right, respectively. In the both cases, we can estimate 
the time delay and Doppler frequency using these contour plots. The proper time delay 
and Doppler difference can be detected by locating the point of symmetry of the 
boomerang. The Figure 37 shows a contour plot when is —0.02 Hz. Since a, is very 
small, the contour plot looks like a straight line. In this case, the dilTerential time delay 
and Doppler difierence can still be estimated by locating the point of symmetry of the 
boomerang ( i.e., center point ). Figure 38 shows the contour plot, when channel- 1 
contains noise only and channel-2 contains signal plus noise ( i.e.. SXR = 10 dB ). Fig- 
ure 39 shows the contour plot, when channel-2 contains noise only and channel- 1 con- 
tains signal and noise ( i.e.. SSR = 10 dB ). Figure 40 shows the contour plot, when 
both channels contain noise only. Using the modified cross power spectrum technique, 
the three types of noise only set ups appear as similiar plots ( see Figures 38. 39 and 
40 ). Figure 38 seems to indicate the presence of signals in both channels, while Figures 
39 and 40 appear more like noise. Hence, contour plots of the modified cross power 
spectrum do not allow clear identification. In this case, it will be difficult to distinguish 
between the noise only situation and a small o <2 situation. 

To use the modified cross power spectrum algorithm at low SXRs ( i.e.. 
SXR < 70 ). we need to investigate the contour plot some more. 
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SNR=10 DEU\Y=0 a1=a2=.23 



o 




Fisjure 35. Contour plot of the modified cross poner spectrum (case 1). 



DOPPLER AXIS 



SNR=10 DEU^Y=0 a1=«2 




SI XV AVnBQ 



Figure 36. Contour plot of the modified cross power spectrum (case 2). 



DOPPLER AXIS 



SNR=10 a1=.23 a2=-.02 DELAY= 




SIXV AV13a 



Figure 37. Contour plot of the modified cross poner spectrum (case 3). 
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DOPPLER AXIS 



SNR=10 CHANNEL1=N0ISE 




Figure 38. Contour plot of the modified cross power spectrum (noise in cliaunel- 1 
only). 
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DOPPLER AXIS 



SNR=10 CHANNEL2=N0ISE 




Figure 39. Contour plot of tlie modified cross power spectrum (noise in cliannel-2 
only). 



66 



DOPPLER AXIS 



CHANNEL2=N0ISE CHANNEL1 =NOISE 



O 




Figure 40. Contour plot of tlie modified cross po>\er spectrum (noise only). 
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DOPPLER AXIS 



APPENDIX E. MODIFIED CROSS POWER SPECTRUM PROGRAM 



Q Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr tV Vr V: Vr :V /V Vr Vr tV Vr */r Vr Vr Vr iV Vr * Vr Vr * Vr Vr Vr Vr : 

C 

main program of approach one. 

this program share subroutine with 
the main program approach two 



* 

* 

* 

Vr 



COMPLEX X(-128; 255), Y( -128: 255), A(0; 2500), B(0: 2500) 
COMPLEX BPFK-128: 191) ,BPF2( -128: 191),BPF(0: 100) 
INTEGER SNR, BINDEL, DELAY 

DATA X,Y,A,B/5770*(0. , 0. ) / ,BPF1 , BPF2 , BPF/741*( 0. ,0. )/ 
OPEN(UNIT=7,FILE=' FILENAME FILETYPE A') 

PI=ACOS(-l. ) 



input part 


f s 


- sampling frequency 


f 


- carrier frequency 


alphal 


- doppler effect at sensorl 


alpha2 


- doppler effect at sensor2 


delay 


- signal receiving time difference 




between two sensors 


SNR 


- signal to noise ratio 


iseed 


- noise generator seed(odd number) 


n 


- number of data during one second 


FS=64. 




F=23. 




ALPHA1=. 


23 


ALPHA2=- 


. 02 


DELAY=0 




SNR=100 




N=64 




T=2. ’’^PI---- 


F/FS 


ISEED=13 




input 


data sampling at two sensors 



CALL SI GNAL( X , F+ALPHA 1 , FS , 0 , SNR , I SEED , 0 , 0 ) 

CALL S IGNAL( Y , F+ALPHA2 , FS , DELAY , SNR , I SEED ,0,0) 
DO 100 K=-128,191 
DO 50 J=K,K+63 



6S 



c * 

c * 

c linear transformation using band pass filter * 

c * 

c * 

OMEGA=T*FLOAT(J-K) 

BPF1(K)=BPF1(K)+X(J)*CMPLX(C0S(0MEGA),-SIN( OMEGA)) 

50 BPF2(K)=BPF2(K)+Y(J)*CMPLX(C0S(0MEGA),SIN(0MEGA)) 

c * 

c * 

c normalization of BPF2 * 

c * 

c - - * 

100 BPF2(K)=BPF2(K)/CABS(BPF2(K)) 

DO 300 BINDEL=-128,127 
DO 200 K=0,N 

c * 

c " 

c input of AR model * 

c * 

c •' 

200 BPF(K)=BPF1(K)*BPF2(K+BINDEL) 

CALL CLEAR(A,B,2500) 

CALL BURGARC BPF ,N+1 , A , IP , VAR) 

CALL CHANGEFFT(B,A,11,MAX) 

DO 250 I=-32,-l 

250 V;RITE(7,997)BINDEL,FL0AT(I)/32. ,CABS(1. /B( 2048+1 ) )*VART 

*,CABS(1. /B(2048+D) 

DO 300 1=0,32 

300 WRITE(7,997)BINDEL,FLOAT(I)/32. ,CABS( 1. /B( I) )*VART 

*,CABS(1. /B(D) 

CL0SE(7) 

STOP 

997 F0RMAT(1X,I4,2X,E13. 6,2X,E13. 6,2X,E13. 6) 

END 
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APPENDIX F, COHERENCE PROGRAM 



c' 

c 

c 

c 

C 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



Vc^V Vv V? Vc Vc^Wc Vc 



main program of approach two 



COMPLEX X(-128: 255) ,Y(-128: 255) ,A(0: 2500) ,B(0: 2500) 
COMPLEX BPFK-128: 191) ,BPF2( -128: 191) ,BPF(0: 100) 
COMPLEX AUT0X(0: 100) , AUT0Y( 0: 100),SUMXY 
INTEGER SNR, BINDEL, DELAY 

DATA X,Y,A,B/5770*(0. ,0. ) / , BPFl , BPF2 , BPF/74l*( 0. ,0. )/ 
DATA AUT0X,AUT0Y/202^’^(0. ,0. )/ 

0PEN(UNTT=7,FILE=’ FILENAME FILETYPE A') 

PI=AC0S(-1. ) 



input part 



f s 
f 

alphal 

alpha2 

delay 

SNR 

iseed 

n 



sampling frequency 
carrier frequency 
doppler effect at sensorl 
doppler effect at sensor2 
signal receiving time difference 
between two sensors 
signal to noise ratio 
noise generator seed(odd number) 
number of data during one second 



ALPHA1=. 23 
ALPHA2=-. 02 
SNR=20 
F=23. 

N=64 

FS=64. 

T=2. ‘-'PI-^T/FS 

ISEED=23 

DELAY=0 



c 

c 

c 

c 

c 



c 

c 

c 

c 

c 



input data sampling at two sensors 



CALL S I GNAL( X , F+ALPHA 1 , FS , 0 , SNR , I SEED , 0 , 0 ) 

CALL S IGNAL( Y , F+ALPHA2 , FS , DELAY , SNR , I SEED ,1,0) 
DO 100 K=-128,191 
DO 50 J=K,K+63 



linear transformation using band pass filter 



v? 

* 



Vf 
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OMEGA=T^'^FLOAT(J-K) 

BPF1(K)=BPF1(K)+X(J)*CMPLX(C0S(0MEGA) ,-SIN(OMEGA)) 
50 BPF2(K)=BPF2(K)+Y(J)*CMPLX(COS(OMEGA) ,SIN(OMEGA)) 

100 CONTINUE 
SUMX=0. 

DO 150 K=0,N 



c 

c find AR order at sensorl * 

c * 

c * 

150 AUT0X(K)=BPF1(K) 

CALL CLEAR(A,B,2500) 

CALL BURGARC AUTOX , N+ 1 , A , I P , VARX ) 

CALL CHANGEFFTC B , A , 1 1 , MAX ) 

XMAX=CABS(1. /B(MAX)) 

DO 300 BINDEL=-128,127 
DO 200 K=0,N 

c -V 

c * 

c input of AR model * 

c * 

c * 

AUTOY(K)=BPF2(K+BINDEL) 

200 BPF(K)=BPF1(K)*(BPF2(K+BINDEL)) 

c * 

c 

c find auto power spectrum of sensor2 at delay k * 

c * 

C -.v 

CALL CLEAR(A.B,2500) 

CALL BURG(AUT0Y,N+1,A,IP,VARY) 

CALL CHANGEFFT(B,A,11,MAX) 

YMAX=1. /CABS(B(MAX)) 

c * 



c 

c find auto power spectrum of cross term •• 

c * 



CALL CLEAR(A,B,2500) 

CALL BURG(BPF,N+1,A,IP,VARXY) 

CALL CHANGEFFTC B, A, 11, MAX) 
XSMAX=(XMAX^>YMAX)*’-'2*VARX*VARY/VARXY*fs 
DO 250 I=-32,-l 

250 WRITE(7,997)BINDEL,FLOAT(I)/32. ,CABS(1. /B( 2048+1 )) /SQRT(XSMAX)*8. 
*,CABS(1. /B(2048+I)) 

DO 300 1=0,32 

300 WRITE( 7,997)BINDEL,FLOAT(I)/32. ,CABS(1. /B( I) )/SQRT(XSMAX)*8. 

*,CABS(1. /B(I)) 

CL0SEC7) 

STOP 



997 F0RMAT(1X,I4,2X,E13. 6,2X,E13. 6,2X,E13. 6) 

END 
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* 

* 



SIGNAL GENERATER * 

INPUT * 

N - NUMBER OF DATA POINTS * 

F1,F2 - FREQUENCY OF 1ST, 2ND SIGNAL 

AMP 1,AMP2 -FREQUENCY OF 1ST, 2ND SIGNAL 

OUTPUT 

A(N) - SIGANL 



SUBROUTINE S IGNAL( B , F , FS , D , SNR , I SEED , ID , NOI SE ) 
COMPLEX B( -128: 255) 

REAL A(0: 1) 

INTEGER SNR,D 
PI=ACOS(-l. ) 

DATA A/0. ,0. / 

SIGMA=1. 

A(0)=SQRT(2.*10. **(SNR/10)) 

IF(NOISE. EQ. 1) A(0)=0. 

DO 100 I=-128,255 
M=(I-D)/128 

IF(I.LT. D. OR.M. NE. 0) M=1 
T=AMOD(F-'-FLOAT( I-D) ,FS) 

CALL GAUSSdSEED, SIGMA, 0. , RANDOM) 

IF(ID.EQ. 0) THEN 

X=C0S( 2. *PI*T/FS)*A(M)+RANDOM 

ELSE 

X=SIN( 2. *PI*T/FS)*A(M)+RANDOM 
ENDIF 

B(I)=CMPLX(X,0. ) 

RETURN 

END 



BURG ALGORITHM 








INPUT 






* 




N 


NUMBER OF DATA POINTS 


* 




X 


INPUT SIGNAL 


* 


OUTPUT 






* 




IP 


ORDER OF AR 


* 




A(0: IP)- 


AR COEFFICIENTS 


•jV 




VAR 


DRIVING NOISE VARIATION 


■>v 






SUBROUTINE BURGAR(X,N,A, IP,VAR) 

COMPLEX X(N) ,A(0: N) ,EFK(500) ,EBK(500) 

COMPLEX EFKK500) ,EBK1(500) ,AA(20,20) ,SUMN,SUTID 
REAL RH0(0: 80) ,FPE(0; 80) 

INTEGER START 
RH0(0)=0. 

FPE(0)=1. E64 
A(0)=CMPLX(1. ,0. ) 

IP=1 
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START=1 
DO 10 1=1, N 

10 RHO(0)=RH0(0)+CABS(X(I))**2/FLOAT(N) 

DO 20 1=2, N 
EFK1(I)=X(I) 

20 EBK1(I-1)=X(I-1) 

LOOP 

K=IP 

SUMN=CMPLX(0. ,0. ) 

SUMD=CMPLX(0. ,0. ) 

DO 30 I=K+1,N 

SUiMN=S UMN+EFK 1 ( I ) *CON JG ( E BK 1 ( I - 1 ) ) 

30 SUMD=SUMD+CABS( EFK1( I )**2 )+CABS( EBK1( I - 1 )**2 ) 

CO SUMD=SUMD+CABS(EFK1( I) )**2+CABS(EBKl( I-l) )**2 

AA(K,K)=-2. *SUMN/SUMD 
TEMP=1. -CABS(AA(K,K))**2 

IF(TEMP. LE. 0. ) TEMP=1. E-10 
RHO ( K ) =TE MP*RHO ( K - 1 ) 

IF(K. GT. 1) THEN 
DO 40 J=1,K-1 

40 AA(J,K)=AA(J,K-1)+AA(K,K)*C0NJG(AA(K-J,K-1)) 

ENDIF 

DO 60 I=K+2,N 

EFK(I)=EFK1(I)+AA(K,K)*EBK1(I-1) 

60 EBK(I-l)=EBKl(I-2)+C0NJG(AA(K,K))*EFKl(I-l) 

DO 70 I=K+2,K 
EFK1(I)=EFK(I) 

70 EBK1(I-1)=EBK(I-1) 

IF(N-K. EQ. 1) THEN 

FPE(K)=FPE(K-1)+1. 

ELSE 

FPE ( K) =RHO( K ) *FLOAT( N+ 1+K ) /FLOAT( N - 1 -K ) 

ENDIF 

IP=IP+1 

UNTIL( FPE(K). GT. FPE(K-l) . AND. K. GT. START) 

IP=K-1 

DO 100 1=1, IP 
100 A(I)=AA(I,IP) 

VAR=RHO(IP) 

RETURN 

END 



C * 

c * 

C BURG ALGORITHM 

C INPUT * 

C N - NUMBER OF DATA POINTS * 

C X - INPUT SIGNAL * 

C IP - ORDER OF AR 

C OUTPUT * 

C A(0: IP)- AR COEFFICIENTS * 

C VAR - DRIVING NOISE VARIATION * 

C * 

C * 
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SUBROUTINE BURG(X,N,A,IP,VAR) 

COMPLEX X(N) ,A(0: N) ,EFK(500) ,EBK(500) 

COMPLEX EFKK500) ,EBK1(500) ,AA(20,20) ,SUMN,SUMD 
REAL RH0(0: 80) ,FPE(0; 80) 

INTEGER START 
RH0(0)=0. 

A(0)=CMPLX(1. ,0. ) 

DO 10 1=1, N 

10 RHO(0)=RHO(0)+CABS(X(I))** 2/FLO AT(N) 

DO 20 1=2, N 
EFK1(I)=X(I) 

20 EBK1(I-1)=X(I-1) 

DO 1000 K=1,IP 

SUMN=CMPLX(0. ,0. ) 

SUMD=CMPLX(0. ,0. ) 

DO 30 I=K+1,N 

SUMN=SUMN+EFK1( I )*CONJG( EBK1( I -1 ) ) 

30 SUMD=SUMD+CABS( EFKK I ) )**2+CABS( EBK1( I - 1) )**2 

AA(K,K)=-2. ’•^SUMN/SUMD 
TEMP=1. -CABS(AA(K,K) )’-*2 

IF(TEMP. LE. 0. ) TEMP=1. E-10 
RH0(K)=TEMP-->RH0(K-1) 

IF(K. GT. 1) THEN- 
DO 40 J=1,K-1 

40 AA( J,K)=AA( J,K-l)+AA(K,K)*CONJG(AA(K-J,K-l)) 

END IF 

DO 60 I=K+2,N 

EFK(I)=EFK1(I)+AA(K,K)*EBK1(I-1) 

60 EBK(I-l)=EBKl(I-2)+CONJG(AA(K,K))*EFKl(I-l) 

DO 70 I=K+2,N 
EFK1(I)=EFK(I) 

70 EBK1(I-1)=EBK(I-1) 

1000 CONTINUE 

DO 100 1=1, IP 
100 A(I)=AA(I,IP) 

VAR=RHO(IP) 

RETURN 

END 



CHANGE FFT 

INPUT 






Vc 




B 


- AR COEFFICIENTS 


* 




A 


- TEMPORARY USING ARRAY 


Vf 


OUTPUT 


ISIZE 


- ALOGCFFT DATA POINTS) /ALOG2 


* 




A 


- POWER SPECTRUM 


Vf 

V? 
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SUBROUTINE CHANGEFFT( A , B , I S I ZE , K ) 
COMPLEX A(0: 2’>*ISIZE- 1) ,B(0: 2**ISIZE-1) 
INTEGER REVRSE 
PI=AC0S(-1. ) 

N=2**ISIZE-1 

K=0 

FS=FL0AT(N+1) 

CALL FFT(N,ISIZE,A,B,FS) 

CALL FINDMAX(A,N,K) 

RETURN 

END 



* 



FINDMAX 






* 


INPUT 






* 




N 


- NUMBER OF SPECTRUM POINTS (N+1) 


Vr 




A 


- POWER SPECTRUM 


* 


OUTPUT 






■jV 




MAX 


- ARRAY INDEX OF MAXIMUM POWER SPECTRUM 


Vf 



SUBROUTINE FINDMAX( A,N,MAX) 
COMPLEX A(0:N) 

MAX=0 

AMAX=1.E66 



DO 100 1=0, N 

IF(CABS(A(D). LT. AMAX) THEN 
MAX=I 

AMAX=CABS(A(I)) 

ENDIF 

CONTINUE 

RETURN 

END 



* 



REVERSE 

BIT REVERSE ORDER INDEX CHANGING SUBROUTINE 
INPUT 

N - ARRAY INDEX 

ISIZE - ALOGCFFT DATA POINTS)/ALOG2 

OUTPUT 

REVERSE- BIT REVERSE ORDER INDEX 



INTEGER FUNCTION REVRSE(N, ISIZE) 
INTEGER A(20) 

DO 200 1=1, ISIZE 
A(I)=MOD(N,2) 

N=N/2 
200 CONTINUE 

REVRSE=0 

DO 300 1=1, ISIZE 
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300 



200 

300 

400 

500 



REVRSE=REVRSE+A( I )*2**( ISIZE-I) 
CONTINUE 
RETURN 
END 



DFT * 

COOLEY - TUKEY ALGORITHM * 

INPUT * 

N - NUMBER OF DATA POINTS * 

ISIZE - AL0G(N)/AL0G2 

A - BIT REVERSE ORDER INDEXED TIME DOMAIN * 

B - TEMPORARY ARRAY 

OUTPUT * 

A - FREQUENCY DOMAIN * 



SUBROUTINE DFT(N , ISIZE , A , B) 

COMPLEX A(0:N),B(0:N),W 
PI=AC0S(-1. ) 

DO 500 11=1, ISIZE 
ISAGE1=2*^HI1-1) 

ISTAGE=2**I1 
DO 200 1=0, N 

ITEST=MOD(I,ISTAGE) 

IFCITEST. GT. ISAGEl) THEN 
K=ITEST-ISAGE1 

THEATA=2. *PI*FLOAT(K)/FLOAT( ISTAGE) 
T1=SIN(THEATA) 

T2=C0S(THEATA) 

IF(ABS(T1). LT. lE-5) T1=0 
IF(ABS(T2). LT. lE-5) T2=0 
V=CMPLX(T2,-T1) 

A(I)=A(I)*W 

END IF 
CONTINUE 
DO 300 1=0, N 

ITEST=MOD(I,ISTAGE) 

IFCITEST. GE. ISAGEl) THEN 

B(I)=-A(I)+A(I-ISAGE1) 

ELSE 

B(I)=A(I)+A(I+ISAGE1) 

ENT) IF 
CONTINUE 
DO 400 1=0, N 
A(I)=B(I) 

CONTINUE 

CO.NTINUE 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



FFT 

1. INDEX CHANGE 

2. USE COOLEY TUKEY ALGORITHM 



INPUT 




N 


- NUMBER OF DATA POINTS 




ISIZE 


- ALOG(N)/ALOG2 




B 


- INPUT(TIME DOMAIN) 




FS 


- SAMPLING FREQUENCY 


OUTPUT 




A 


- FREQUENCY DOMAIN 



* 

* 



* 



SUBROUTINE FFT(N,ISIZE,A,B,FS) 
COMPLEX A(0: N),B(0:N) 

INTEGER REVRSE 
DO 100 1=0, N 

K=I 

100 A(REVRSE(K,ISIZE))=B(I) 

CALL DFT(N,ISIZE,A,B) 

DO 200 1=0, N 

200 A(I)=A(I)/FS 

RETURN 
END 

C 

c 

C CLEAR 



C 


INPUT 








C 




N 


- NUMBER OF DATA 


POINTS 


C 




A,B 


- INPUT arrays 




C 


OUTPUT 








C 




A,B 


- reinitiallized 


arrays 



C 



C 

SUBROUTINE CLEAR(A,B,N) 

COMPLEX A(0;N),B(0;N) 

DO 100 1=0, N 

100 A(I)=B(I)=CMPLX(0. ,0. ) 

RETURN 
END 
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